
# Replication file for Table 2 (body) and Appendix G.2. Tables 18,20,22,24,26.


#install.packages("Hmisc")
require(Hmisc)
#install.packages("stargazer")
require(stargazer)
require(multiwayvcov)
require(lmtest)

load("Brexit_Rep.RData")

treated.sub <-  brexit.sub2[which(brexit.sub2$post == 1),]
control.sub <-  brexit.sub2[which(brexit.sub2$post == 0),]


varlabels <- c("Treat", "Male", "Age", "Income (Imp.)", "Social Grade", "Education","Employed", "Children (=1)", "Married", "North West", 
               "Yorkshire and the Humber", "East Midlands", "West Midlands", "East of England","London", "South East", 
               "South West", "Wales", "Scotland", "Constant")

controls.names <- Cs(Gender, Age, Hincome.imp, socialgrade, Education, employed, children, marstat, GOR_full)

right.side.eq <- "~ wave + Gender + Age + Hincome.imp + socialgrade + Education + employed + children + marstat + as.factor(GOR_full)"


### Control Wave 6 and Control Wave 7 (June 2016, Pre Referendum/ Nov. 2016)


ones <- rep_len(1, length.out = nrow(control.sub))
twos <- rep_len(2, length.out = nrow(control.sub))


dvs.6 <- cbind(ones, control.sub$w8w6, control.sub$leaver, control.sub$migs.take.jobs.6, control.sub$migs.more.terror.6,control.sub$closed.immigrants.6,
               control.sub$hurt.standing.refugee.6,
               control.sub$threaten.culture.refugee.6,control.sub$overwhelm.welfare.refugee.6, control.sub$GUNQID)

dvs.7 <- cbind(twos,control.sub$w8w6, control.sub$leaver, control.sub$migs.take.jobs.7,control.sub$migs.more.terror.7,control.sub$closed.immigrants.7,
               control.sub$hurt.standing.refugee.7,
               control.sub$threaten.culture.refugee.7,control.sub$overwhelm.welfare.refugee.7, control.sub$GUNQID)

controls <- cbind(control.sub$Gender, control.sub$Age, control.sub$Hincome.imp, control.sub$socialgrade, control.sub$Education, 
                  control.sub$employed, control.sub$children, control.sub$marstat, control.sub$GOR_full)

dvs <- Cs(migs.take.jobs,migs.more.terror,closed.immigrants,hurt.standing.refugee,
                threaten.culture.refugee,overwhelm.welfare.refugee)

dvs.6 <- cbind(dvs.6, controls)
dvs.7 <- cbind(dvs.7, controls)

c67 <- as.data.frame(rbind(dvs.6, dvs.7))


names(c67) <- c("wave", "w8w6", "leaver",dvs,"GUNQID",controls.names)

c67.l <- subset(c67, leaver == 1)
c67.r <- subset(c67, leaver == 0)

lf.l <- vector(length(dvs), mode = "list")
names(lf.l) <- dvs

lf.r <- vector(length(dvs), mode = "list")
names(lf.r) <- dvs


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  x <- eval(substitute(lm(.modelformula, data = c67.l, 
                                        weights = w8w6), list(.modelformula = modelformula)))
  

  
 x.clust <- cluster.vcov(x, c67.l$GUNQID)
 x.mcse <- coeftest(x, x.clust)
  
  
lf.l[[dvs[i]]] <- x.mcse
  
print(lf.l[[dvs[i]]])

}


for (i in 1:(length(dvs))){
  
  modelformula <- paste(dvs[i],right.side.eq)
  print(modelformula)
  
  x.r <- eval(substitute(lm(.modelformula, data = c67.r, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  x.clust.r <- cluster.vcov(x.r, c67.r$GUNQID)
  x.mcse.r <- coeftest(x.r, x.clust.r)
  
  lf.r[[dvs[i]]] <- x.mcse.r
  
  print(lf.r[[dvs[i]]])
  
}


# Refugees

stargazer(lf.l[[dvs[5]]], lf.r[[dvs[5]]], lf.l[[dvs[6]]], lf.r[[dvs[6]]],lf.l[[dvs[4]]], lf.r[[dvs[4]]], covariate.labels = varlabels)

# Migrants

stargazer(lf.l[[dvs[3]]], lf.r[[dvs[3]]], lf.l[[dvs[1]]], lf.r[[dvs[1]]],lf.l[[dvs[2]]], lf.r[[dvs[2]]], covariate.labels = varlabels)


#### Treated Wave 5 and Wave 6  (Dec. 2015/ July 2016, Post Referendum)

ones <- rep_len(1, length.out = nrow(treated.sub))
twos <- rep_len(2, length.out = nrow(treated.sub))

t.ref.test.5 <- cbind(ones,treated.sub$w8w6,treated.sub$leaver,treated.sub$hurt.standing.refugee.5,
                      treated.sub$threaten.culture.refugee.5,treated.sub$overwhelm.welfare.refugee.5, 
                      treated.sub$closed.immigrants.5, treated.sub$GUNQID)
t.ref.test.6 <- cbind(twos,treated.sub$w8w6,treated.sub$leaver,treated.sub$hurt.standing.refugee.6,
                      treated.sub$threaten.culture.refugee.6,treated.sub$overwhelm.welfare.refugee.6,
                      treated.sub$closed.immigrants.6, treated.sub$GUNQID)

controls <- cbind(treated.sub$Gender, treated.sub$Age, treated.sub$Hincome.imp, treated.sub$socialgrade, treated.sub$Education, 
                  treated.sub$work_stat_full, treated.sub$children, treated.sub$marstat, treated.sub$GOR_full)

t.ref.test.5 <- cbind(t.ref.test.5, controls)
t.ref.test.6 <- cbind(t.ref.test.6, controls)

t.w56 <- as.data.frame(rbind(t.ref.test.5, t.ref.test.6))

ref.names <- c("hurt.standing.refugee","threaten.culture.refugee","overwhelm.welfare.refugee","closed.immigrants")

names(t.w56) <- c("wave", "w8w6", "leaver",ref.names, "GUNQID",controls.names)

t.w56.l <- subset(t.w56, leaver == 1)
t.w56.r <- subset(t.w56, leaver == 0)

tl <- vector(length(ref.names), mode = "list")
names(tl) <- ref.names

tr <- vector(length(ref.names), mode = "list")
names(tr) <- ref.names

#Leavers
for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
  x.l <- eval(substitute(lm(.modelformula, data = t.w56.l, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  
  
  x.clust.l <- cluster.vcov(x.l, t.w56.l$GUNQID)
  x.mcse.l <- coeftest(x.l, x.clust.l)
  
  
  tl[[ref.names[i]]] <- x.mcse.l
  
  print(tl[[ref.names[i]]])
  
}

#Remainers

for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
 x.r <- eval(substitute(lm(.modelformula, data = t.w56.r, 
                                       weights = w8w6), list(.modelformula = modelformula)))
  
  
 x.clust.r <- cluster.vcov(x.r, t.w56.r$GUNQID)
 x.mcse.r <- coeftest(x.r, x.clust.r)

 
 
  tr[[ref.names[i]]] <- x.mcse.r
  print(tr[[ref.names[i]]])
  
}


stargazer(tl[[ref.names[2]]], tr[[ref.names[2]]], tl[[ref.names[3]]], tr[[ref.names[3]]], tl[[ref.names[1]]], tr[[ref.names[1]]], tl[[ref.names[4]]], tr[[ref.names[4]]], covariate.labels = varlabels)



##Control Wave 6 and Control Wave 8 (Pre-Referendum: June 2016/ August 2017)

ones <- rep_len(1, length.out = nrow(control.sub))
twos <- rep_len(2, length.out = nrow(control.sub))


dvs <- Cs(closed.immigrants,hurt.standing.refugee, threaten.culture.refugee,overwhelm.welfare.refugee)

dvs.6 <- cbind(ones, control.sub$w8w6, control.sub$leaver, control.sub$closed.immigrants.6,
               control.sub$hurt.standing.refugee.6, control.sub$threaten.culture.refugee.6,
               control.sub$overwhelm.welfare.refugee.6,control.sub$GUNQID)

dvs.8 <- cbind(twos, control.sub$w8w6, control.sub$leaver, control.sub$closed.immigrants.8,
               control.sub$hurt.standing.refugee.8, control.sub$threaten.culture.refugee.8,
               control.sub$overwhelm.welfare.refugee.8,control.sub$GUNQID)


controls <- cbind(control.sub$Gender, control.sub$Age, control.sub$Hincome.imp, control.sub$socialgrade, control.sub$Education, 
                  control.sub$employed, control.sub$children, control.sub$marstat, control.sub$GOR_full)


dvs.6 <- cbind(dvs.6, controls)
dvs.8 <- cbind(dvs.8, controls)


c.w68 <- as.data.frame(rbind(dvs.6, dvs.8))

names(c.w68) <- c("wave", "w8w6", "leaver",ref.names,"GUNQID",controls.names)

c.w68.l <- subset(c.w68, leaver == 1)
c.w68.r <- subset(c.w68, leaver == 0)


cl <- vector(length(ref.names), mode = "list")
names(cl) <- ref.names

cr <- vector(length(ref.names), mode = "list")
names(cr) <- ref.names

#Leavers

for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
 x.l <- eval(substitute(lm(.modelformula, data = c.w68.l, 
                                           weights = w8w6), list(.modelformula = modelformula)))
  
 x.clust.l <- cluster.vcov(x.l, c.w68.l$GUNQID)
 x.mcse.l <- coeftest(x.l, x.clust.l)
 
 
 cl[[ref.names[i]]] <- x.mcse.l
 
 print(cl[[ref.names[i]]])
  
}

#Remainers
for (i in 1:(length(ref.names))){
  
  modelformula <- paste(ref.names[i],right.side.eq)
  print(modelformula)
  
  x.r <- eval(substitute(lm(.modelformula, data = c.w68.r, 
                                           weights = w8w6), list(.modelformula = modelformula)))
  
  
  x.clust.r <- cluster.vcov(x.r, c.w68.r$GUNQID)
  x.mcse.r <- coeftest(x.r, x.clust.r)
  
  cr[[ref.names[i]]] <- x.mcse.r
  
  print(cr[[ref.names[i]]])
  
}


stargazer(cl[[ref.names[2]]], cr[[ref.names[2]]], cl[[ref.names[3]]], cr[[ref.names[3]]], cl[[ref.names[1]]], cr[[ref.names[1]]], cl[[ref.names[4]]], cr[[ref.names[4]]], covariate.labels = varlabels)


